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ABSTRACT 

We introduce a new class of 2-D mass models, whose potentials are of Stackel form 
in elliptic coordinates. Our model galaxies have two separate strong cusps that form 
double nuclei. The potential and surface density distributions are locally axisymmetric 
near the nuclei and become highly non- axisymmetric outside the nucleus. The surface 
density diverges toward the cuspy nuclei with the law E oc r~^. Our model is sustained 
by four general types of regular orbits: butterfly, nucleuphilic banana, horseshoe and 
aligned loop orbits. Horseshoes and nucleuphilic bananas support the existence of 
cuspy regions. Butterflies and aligned loops control the non-axisymmetric shape of 
outer regions. Without any need for central black holes, our distributed mass models 
resemble the nuclei of M31 and NGC4486B. It is also shown that the self-gravity of 
the stellar disc can prevent the double nucleus to collapse. 
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1 INTRODUCTION 

Hubble Space Telescope (HST) data revealed that M31 and 
NGC4486B have double nuclei (Lauer et al. 1996, hereafter 
L96; Tremaine 1995, hereafter T95). M31 has a bright nu- 
cleus (PI) displaced from the centre of the isophotal lines of 
outer regions and a fainter nucleus (P2) just at the centre. 
NGC4486B exhibits a similar structure with a minor dif- 
ference: The centre of outer isophotes falls between PI and 
P2. There are some explanations for the emergence of the 
double nuclei of these galaxies, among which the eccentric 
disc model of T95 has been more impressive. In the model 
of T95, a central black hole (BH) enforces stars to move 
on "aligned" Keplerian orbits, which may elongate in the 
same direction as the long-axis of the model. Stars moving 
on aligned Keplerian orbits linger near apoapsis and may 
result in PI. The mass of central "supermassive" BH should 
be much greater than the mass of neighboring disc. Other- 
wise, the asymmetric growth of PI won't allow the BH to 
remain in equilibrium. 

Goodman & Binney (1984) showed that central mas- 
sive objects enforce the orbital structure of stellar systems 
to evolve towards a steady symmetric state. This result was 
then confirmed by the findings of Merritt & Quinlan (1998) 
and Jalali (1999, hereafter J99) in their study of elliptical 
galaxies with massive nuclear BHs. Within the BH sphere of 
influence, highly non-axisymmetric structure can only exist 
for a narrow range of BH mass (J99). The results of J99 
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show that long-axts tube orbits of non-axisymmetric discs 
with central massive BHs, elongate in the both ± directions 
of long-axis. Thus, the probabilities for the occurrence of 
two bright regions, in both sides of BH along the long-axis, 
are equal (these bright regions are supposed to be formed 
near the apogee of long-axis tubes) . By this hypothesis one 
can interpret the double structure of NGC4486B by placing 
a supermassive BH between PI and P2. However, some dis- 
advantages arise in the case of M31. In the nucleus of M31, 
the formation of PI can still be deduced from the behav- 
ior of long-axis tubes. But, there is no mathematical proof 
for the "coexistence" of PI and P2 when the centre of P2 
coincides with BH's location. 

In this paper we attempt to create a model based on 
the self-gravity of stellar discs to show that systems with 
double nuclei can exist even in the absence of central BHs. 
Especial cases of our non-scale-free planar models are eccen- 
tric discs, which display a collection of properties expected 
in self-consistent cuspy systems. Our models are of Stackel 
form in elliptic coordinates (e.g., Binney & Tremaine 1987) 
for which the Hamilton- Jacobi equation separates and stel- 
lar orbits are regular. 

In most galaxies, density diverges toward the centre in 
a power-law cusp. In the presence of a cusp, regular box 
orbits are destroyed and replaced by chaotic orbits (Ger- 
hard & Binney 1985). Through a fast mixing phenomenon, 
stochastic orbits cause the orbital structure to become ax- 
isymmetric at least near the centre (Merritt & Valluri 1996). 
These results are confirmed by the findings of Zhao et al. 
(1999, hereafter Z99). Their study reveals that highly non- 
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axisymmetric, scale- free mass models can not be constructed 
self-consistently. Near the cuspy nuclei, the potential func- 
tions of our distributed mass models are proportional to 
as r ^ 0. So, we attain an axisymmetric structure near the 
nuclei which is consistent with the mentioned nature of den- 
sity cusps. The slope of potential function changes sign as 
we depart from the centre and our model galaxies consid- 
erably become non-axisymmetric. Non-axisymmetric struc- 
ture is supported by butterfly and aligned loop orbits. Close 
to the larger nucleus, loop orbits break down and give birth 
to a new family of orbits, horseshoe orbits, which in turn gen- 
erate nucleuphilic banana orbits. Stars moving in horseshoe 
and banana orbits lose their kinetic energy as they approach 
to the nuclei and contribute a large amount of mass to form 
cuspy regions. 



2 THE MODEL 

Let us introduce the following family of planar potentials 
expressed in the usual [x, y) cartesian coordinates: 

-1 



(ri + r2)^ - /9(ri - r2)\ri - r2p 



2rir2 

^-{x + af + y"^, rl = {x- 



n2 , 2 

a) +y , 



(1) 
(2) 



where a, K>0, 0<[3<1 and 2 < 7 < 3 are constant 
parameters. The points (a; = —a,y = 0) and {x = a,y = 0) 
are the nuclei of our 2-D model. We call them PI and P2, 
respectively. The distance between PI and P2 is equal to 
2a. The surface density distribution corresponding to $ is 
determined by (Binney & Tremaine 1987) 

(V^$)da;d?/ 



E{x',y') 



1 



(3) 



It is a difficult task to evaluate (|3|) analytically. So, we have 
adopted a numerical technique to calculate this double in- 
tegral. The functions $ and E are cuspy at PI and P2. To 
verify this, we investigate the behavior of $ and E near the 
nuclei (ri and r2 0). Sufficiently close to PI, we have 
ri <^ r2 that simplifies (nl) as follows 



KrT^il+^r+fil^ 
2 ri 

We expand (1 + and (1 - \ 

obtain 



(4) 

in terms of ri/r2 to 



$ : 



Kr: 



7-1 



2ri 



(l + (-l)"/3)r(7 + l) fn 



n!r(7 -n+1) 



(5) 



where F is the well known Gamma function. As n tends to 
zero, r2 is approximated by 2a and ri/r2 0. Therefore, 
Equation (^) reads 

0(1 



Q = -K(2a 



(6) 



ri 2 
from which one concludes 

Ecx(3(l+/3)r-i"^ (7) 

Similarly, it can readily be shown that the following approx- 
imations hold close to P2 (tiItx — > 0), 

0(1-/3) 



ri 



(8) 



E cx Q(l-/3)r2-'. (9) 

In distant regions, when r ^ a (with = -|- y^), the 
potential function is approximated by 

$ K. 2''-^Kr"'-^. 

Correspondingly, 



(10) 



E oc r 



,7-3 



(11) 



This shows that the surface density falls off outward if 7 < 3. 
Besides, orbits will be bounded if the potential $ is concave 
in outer regions. This requirement implies 7 > 2. Thus, we 
are restricted to 2 < 7 < 3. In Fig. 1, we have plotted the 
isocontours of $ and E for 7 = 2.8, a — 0.5 and several 
choices of l3. The 3-D views of $ and E have also been 
demonstrated in Fig. 2. In § 5, the potential surface of Fig. 2a 
will be referred as potential hill. 

Figs. 1 and 2 assure that the potential and surface den- 
sity functions are cuspy at PI and P2. Regardless of the 
values of constant parameters, the potential $ has a local 
minimum at {x = 0,y = 0). This minimum point can eas- 
ily be distinguished in Figs, la, lb and Ic. As the surface 
density isocontours show, the cuspy zones are disjointed by 
two separatrices that transversally intersect each other at 
a saddle point located on the a;-axis between PI and P2. 
The a;-coordinate of this point can be determined through 
solving 



= 0, y = 0, 



(12) 



d^ix,y) 

dx 
for X. 

The parameter /3 controls the sizes of cuspy zones 
around PI and P2. For (3 = 0, the sizes of two cuspy zones 
are equal and the model has reflection symmetries with re- 
spect to coordinate axes. For < /? < 1, the size of cuspy 
zone near PI is larger than that of P2 and the model is 
only symmetric with respect to the a;-axis. The cuspy region 
around P2 is shrunk to zero size when /3 = 1 and we at- 
tain an eccentric disc with a single nuclear cusp. Equations 
(^) through (P) show that the potential and surface den- 
sity functions are approximately axisymmetric in the neigh- 
borhood of PI and P2. As we move outward, a "highly" 
non-axisymmetric structure occurs. For the large values of 
r, the surface density monotonically decreases outward and 
our model galaxies become rounder again. Our mass models 
are indeed hybrid ones, which reflect the properties of den- 
sity cusps and non-axisymmetric systems, simultaneously. 
The centre of outer surface density isocontours falls at the 
middle of the centerline of PI and P2. Nevertheless, the ef- 
fective cuspy zones around PI and P2 have difi'erent sizes. 

In what follows, we show that the potential $ is of 
Stackel form in elliptic coordinates. We then classify pos- 
sible orbit families, all of which are non-chaotic. 



3 ORBIT FAMILIES 

We carry out a transformation to elliptic coordinates as fol- 
lows 



X = a cosh M cos 
y = asinhttsint). 



(13) 
(14) 
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id) (e) {/) 

Figure 1. The potential (#) and surface density (S) isocontours for 7 = 2.8, a = 0.5 and K = 0.2. Figures (a), (b) and (c) show the 
potential isocontours for /3 = 0, /3 = 0.75 and (3 = 1, respectively. The corresponding surface density isocontours are illustrated in Figures 
(d),(e) and (f). 




(a) (6) 

Figure 2. The 3-D views of $ and S for 7 = 2.8, a = 0.5, K = 0.2 and /3 = 0.75. (a) the potential function (b) the surface density 
distribution. 
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where m > and < « < 27r. The curves of constant it and 
V are confocal ellipses and hyperbolas, respectively. PI and 
P2 are the foci of these curves. In the new coordinates, the 
motion of a test star is determined by the Hamiltonian 



H 



1 



-{Pu +Pv) +"1>(U,U), 



2a2(sinh^ u + sin'^ v) 
where pu and Pv denote the canonical momenta and 
^ _ F{u) + G{v) 



(15) 



(16) 

(17) 
(18) 



(sinh^ u + sin^ v) ' 
F{u) = C(cosh?i)^, 
G(v) = — C/3cosu| cos^P""'", 
C = K{2ay. 

The transformed potential ( |l^ ) is of Stackel form for which 
the Hamilton- Jacobi equation separates and yields the sec- 
ond integral of motion, l2- We obtain 



h=pl- 20^ E smh^ u-^ F{u), 
or equivalently 

- h^pl- 2a^Esm v + G{v) 



(19) 



(20) 



where E is the total energy of the system, E = Ti.. The 
potential function ($) is positive everywhere. Hence, we im- 
mediately conclude E > 0. 



Having the two isolating integrals E and I2, one can find 
the possible regions of motion by employing the positiveness 
of p^ and pi in ( |l^ and (^o|). We define the following func- 
tions: 



f{u) = -2a'' E sinh' u + F{u), 
g{v) = -2a'Esm''v + G{v). 

Since > and > 0, one can write 

h-f{u) > 0, 
~l2 - g{v) > 0. 



(21) 
(22) 

(23) 
(24) 



Orbits are classified based on the behaviour of f{u) and g{v). 
The most general form of f{u) is attained for 7C < 4a'E. 
In such a circumstance, /(u) has a local maximum at « = 0, 
/m = /(O) — C, and a global minimum a,t u = Um, fm ~ 
f{um), where 



cosh Un 



and 



4a'' E \ 



fm = — 2a^iJsinh^ Um + C(coshiim)^. 
From ( p^ we obtain 

I2 > fm. 



(25) 



(26) 



(27) 
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On the other hand, 17(1;) has a global maximum at v = n, 
qm = ~ PC , and two global minima at w = -e 12 

and V = 3n/2, gn^=g{n /2)=g(3n /2)=-2a'^ E. Therefore, In- 
equality (M) implies 



I2 < 2a E. 

By combining (|27|) and one achieves 



(28) 



/m <l2< 2a^E. (29) 

By taking 2 < 7 < 3 and 7C < Aa?E into account, we arrive 
at 2a?'E > C. Furthermore, /m and in consequence I2, can 
take both positive and negative values. For a specified value 
of E, the following types of orbits occur as I2 varies. 



(i) Butterflies. 
for u and v axe 



For C < I2 < 2a E, the aUowed values 



It < Uo, Vb,l < V < Vb,2, Vb,3 < V < Vb,^ 



(30) 



where mo and Vb,i {i — 1,2,3,4) are the roots of f{u) — I2 
and g{v) — —I2, respectively. As Fig. 3a shows, the horizon- 
tal line that indicates the level of I2, intersects the graph 
of f{u) at one point, which specifies the value of wo- The 
line corresponding to the level of —I2 intersects g{v) at four 
points that give the values of Wi,,iS (Fig. 3b). In this case the 
motion takes place in a region bounded by the coordinate 
curves u = uq and v = Vb.i- The orbits fill the shaded region 
of Fig. 4a. These are butterfly orbits (de Zeeuw 1985) that 
appear around the local minimum of "1> at [x = Q,y = G). 

(ii) Nucleuphilic Bananas. For f3C < I2 < C the equa- 
tion f{u) — I2 has two roots, u„,i and u„,2, which can be 
identified by the intersections of f{u) and the level line of I2 
(see Fig. 3c). In this case, the equation g{v) = —I2 has four 
real roots, v = Vn,i {i = 1,2,3,4), (Fig. 3d). The allowed 
ranges of u and v will be 



"U-n,! <U< Un,2, Vn,l <V< V„,2, Vn,3 < V < Vn 



(31) 



The orbits (Fig. 4b) are bound to the curves of u = Mn,i, 
u = M„,2 and V — Vn.i- We call them nucleuphilic banana 
orbits, for they look like banana and bend toward the nuclei. 

(iii) Horseshoes. For —(3C < I2 < PC, both of the equa- 
tions f[u) = I2 and g{v) = —I2 have two real roots. We 
denote these roots by m = Uh,i and v = Vh,i {i = 1,2). In 
other words, the level lines of ±/2 intersect the graphs of 
f{u) and g{v) at two points as shown in Figs. 3e and 3f. 
The trajectories of stars fill the shaded region of Fig. 4c. We 
call these horseshoe orbits. 

(iv) Aligned Loops. For fm < I2 < —PC, the equation 
f{u) = I2 has two real roots, u = ui^t (i = 1,2) while the 
equation g{v) = —I2 has no real roots and Inequality (^|) is 
always satisfied (Figs. 3g and 3h). The orbits fill a tubular 
region as shown in Fig. 4d. We call these aligned loops be- 
cause they are aligned with the surface density isocontours 
of outer regions. 

(v) Transitional cases. For I2 ~ 2a^E, stars undergo 
a rectilinear motion on the j/-axis with the amplitude of 
iasinhiiQ. For I2 ~ fm, loop orbits are squeezed to an el- 
liptical orbit defined by ii = Um. For P — 0, horseshoe orbits 
are absent, leaving the other types of orbits symmetric with 
respect to the coordinate axes. Banana orbits no longer sur- 
vive for P = 1 (eccentric disc model). In this case, butter- 
flies extend to a lens orbit when I2 = C (see Figure 4e). For 
7C > 4a2s, f{u) is a monotonically increasing function of u 
and "low-energy" butterflies are the only existing family of 




(a) 





id) 




Figure 4. The possible families of orbits: (a) a butterfly orbit (b) 
nucleuphilic banana orbits (c) a horseshoe orbit (d) an aligned 
loop orbit (e) a lens orbit associated with (3 = 1 and I2 = C. 
For /9 ^ 0, the orbits are only symmetric with respect to the x- 
axis. Loop orbits are exceptional; they are always symmetric with 
respect to the coordinate axes. 



orbits. These are small-amplitude liberations in the vicinity 
of the local minimum of <1> at (a; = 0, y = 0). 



4 THE POSITIVENESS OF THE SURFACE 
DENSITY 

The sign of E is linked to that of V^^ through Equation (|). 
To prove that E takes positive values for the potentials of 
(0), it suffices to show that the Laplacian of <1> is a positive 
function of 7, /3, u and v. 

Consider the Laplace equation in elliptic coordinates as 



(32) 



D = sinh^ It + sin^ ?;, 



where ,3 denotes Substituting from ( p^ into (|3^, leads 
to 



^ ■^(7,/3;m,w) 



2a*L»* 



(33) 



with 



r = D^{F^uu + G,^^)-D{F + G){D,uu+D,^^) 
-2D{F,uD,u + G,vD,^) 
+2iF + G){D%+D%). 



(34) 
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Figure 5. The behaviour of Q{'y,v) for — < v < ^ and 2 < 
7 < 3. 

For the sake of simplicity, we assume C = 1. We show 
that the minimum of is always positive. We prove our 
claim for — -f < f < f-, which implies G{v) — ~l3cos''v 
(a similar method can be repeated for ^ < v < ^). In 
this case, J- will be a linear, decreasing function of f3 (be- 
cause $ has such a property). Therefore, one concludes 
J^{'y,l;u,v) < J^{'y,j3;u,v). Furthermore, directly de- 
pends on cosh u, which results in ^^(7, 1; 0, v) < J-i-y, 1; u, v). 
Hence, is positive if G{'y, v) = J'i'J, 1; 0, v) > 0. By the 
evaluation of ( |34[ ) for 13 = 1 and u = 0, one finds out 

Q{-f,v) = — 7^ sin^ n cos^"""^ « + sin^ — cos'' u) 
-|-7[sm n -I- cos ii(sm v 
— 3 cos^ u sin* w)]. (35) 

We have plotted 6(7,17) in Figure 5. On the evidence of 
this figure, Q is a positive function for — < w < § and 
2 < 7 < 3. Thus, the surface density distribution takes 
positive values for all of our model galaxies. 



5 DISCUSSION 

In his pioneering work, Euler showed the separability of 
motion in the potential field of two fixed Newtonian cen- 
tres of attraction. This problem was then completely solved 
by Jacobi (Pars 1965). It is physically impossible to keep 
apart these two "point masses", for they will attract each 
other leading to an eventual collapse. However, the assumed 
masses can be in equilibrium if they revolve around their 
common centre of mass (this is the classical 3-body prob- 
lem). Our planar model is indeed Jacobi's problem in which 
we have replaced two fixed centres of gravitation with a con- 
tinuous distribution of matter, where mass concentration in- 
creases towards two nuclei (PI and P2) in power-law, strong 
cusps. These nuclei are maintained by an interesting family 
of orbits, nucleuphilic bananas. Below, we explain why the 
mentioned nuclei are generated and don't collapse. 

The force exerted on a star is equal to — V$. The mo- 
tion under the influence of this force can be tracked on the 
potential hill of Figure 2a. This helps us to better imagine 
the motion trajectories. 



5.1 The behaviour of orbits 

Stars moving in nucleuphilic banana orbits undergo mo- 
tions near the 1:2 resonance. They oscillate twice in the 
{/-direction for each x-axial oscillation. The turning points 
of this group of stars lie on the curves v = Vn,i- These hy- 
perbolic curves can be chosen arbitrarily close to PI or P2. 
When stars approach PI (or P2), their motion slows down 
(because they climb on the cuspy portion of the potential hill 
and considerably lose their kinetic energy while the potential 
energy takes a maximum) and the orbital angular momen- 
tum switches sign somewhere on v = Vn,i- Thus, these stars 
spend much time in the vicinity of PI (or P2) and deposit a 
large amount of mass. This phenomenon is the main reason 
for the generation of cuspy zones around PI and P2. Stars 
moving in nucleuphilic bananas cross the y-axis quickly, and 
therefore, don't contribute much mass to the region between 
the nuclei. 

Horseshoe orbits cause the sizes of cuspy zones to be 
different through the following mechanism. Stars that start 
their motion sufficiently close to PI (larger nucleus), are 
repelled from PI because the force vector is not directed in- 
ward in this region. As they move outward, their orbits are 
bent and cross the x-axis with non-zero angular momen- 
tum. These stars linger only near PI, and in consequence, 
help the cuspy zone around PI grow more than that of P2. 
The asymmetry of nucleuphilic bananas, with respect to the 
y-axis, is also an origin of the different sizes of cuspy zones. 
In fact, horseshoe orbits are born once nucleuphilic bananas 
join together for I2 = PC . Horseshoe and nucleuphilic ba- 
nana orbits are the especial classes of boxlets that appropri- 
ately bend toward the nuclei. The lack of such a property 
in centrophobic banana orbits causes the discs of Sirdhar & 
Touma (1997) to be non-self-consistent. 

Aligned loop orbits occur when the orbital angular mo- 
mentum is high enough to prevent the test particle to slip 
down on the potential hill. The boundaries of loop orbits are 
defined by the ellipses u = and u = «;,2. The nuclear 
cusps are located at the foci of these ellipses. Aligned loops 
have the same orientation as the surface density isocontours 
(compare Figures 1 and 4d). Thus, according to the results 
of Z99, it is possible to construct a self-consistent model 
using aligned loop orbits. 

It is worthy to note that butterfiy orbits play a signifi- 
cant role in maintaining the non-axisymmetric structure of 
the model at the moderate distances of 0{a). 

5.2 The nature of PI and P2 

The points where the cusps have been located, are inherently 
unstable. With a small disturbance, stars located at {x = 
±a,y = 0) are repelled from these points because — V"I> is 
directed outward when ^ (i = 1,2). But, the time 
that stars spend near the nuclei will be much longer than 
that of distant regions when they move in horseshoe and 
banana orbits. The points PI and P2 are unreachable, for 
they correspond to the energy level E = -l-oo. Based on the 
results of this paper, we conjecture that there may not be 
any mass concentration just at the centre of cuspy galaxies. 
However, a very dense region exists arbitrarily close to the 
centre! 
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Figure 6. In this figure we have shown a circular ring of matter 
of outer regions. The masses of cuspy zones can approximately 
be computed and concentrated at the points (±a, 0). The gravi- 
tational attraction of the ring upon the point mass Mi {M2) is 
directed in the x-direction due to the existing symmetry. 



5.3 The double nucleus can be in equilibrium 

The nuclei pull each other due to their mutual gravitational 
attraction and it seems that they must collapse. However, we 
explain that in certain circumstances, the double nucleus can 
be in static equilibrium. At first we estimate the mass inside 
the separatrices of the surface density distribution (the mass 
of cuspy zones) and concentrate the matter at PI and P2 
(this is logical because the surface density distribution is 
almost axisymmetric near the nuclei). In this way, we obtain 
two point masses. Mi and A/2. According to (Q) and (^, the 
following relations approximately hold 



air 



roi 



= 2Tvai log — 



1,2, 



0, 



(36) 



where roi are chosen as the radii of inner tangent circles to 
the separatrices and the constant parameters ai are com- 
puted based on the surface density profile near the nuclei. 
As e — > 0, MiS diverge to infinity unless a negative mech- 
anism prevents them to grow. Consider the discs of radius 
e with the centres located at (±a, 0) and call them Di and 
T>2- Since PI and P2 are "locally" unstable, we rely on our 
previous argument that the matter is swept out from these 
points, allowing us to exclude ©i and ©2 from our model 
for some < e <C 1. In this way, MiS take finite values. 

Ml and M2 attract each other and start to move if they 
are not infiuenced by other gravitational sources. We claim 
that the required extra force comes from the gravitational 
attraction of the matter of outer regions. Consider Figure 6 
where Mi and M2 are shown along with a ring of matter of 
outer regions. For brevity, we assume /3 = 0, which yields 
Ml = M2 = M. Due to the existing symmetry, the gravi- 
tational force exerted on M2 by the assumed ring will have 
a resultant in the a;-direction. When r is sufficiently large, 
r ^ a, this force is calculated as follows 



Mr) 



GMcjo 



.7-2 



cos (f)A6 



2ar cos 6 



(r^+a"^ -2arcos6i)i/2' 



(37) 
(38) 



where we have used E 



with (Too being a positive 



constant (see Eq. (|llD). By integrating Fx{r) over r from 



some r = R^a to r = 00, the total force, due to the 
matter of outer regions, is found to be 



= GMao 



r"' ^ (r cos 9 — a)A9Ar 
(r^ + — 2ar cos 9)^^^ 



(39) 



By a change of independent variable as ^ = a/r, the inte- 
grand can be simplified in the form 

(40) 



Vo ^''^'^^dUo (l + C^-2Ccose)V2 

2GM(Jac 



where (,0 — a/R. Consequently, 



«0 



P„{cos9)d9, 



(41) 



(42) 



with P„s being the well known Legendre functions. Accord- 
ing to (Morse & Feshbach 1953) 



2k + l 



{cos9)d9 = 0, 



P2fc(cos6i)d6i 



(2fe)! 
(2'=fc!)2 



C2fc, 



one achieves 
rio 



/•to 

= 6/ (C'~^^2fcc2fee''"')dC- 
Jo k=l 



Integrating ( |45[ ) over ^, yields 
F. = 6Q(Co), 



Q(eo) = Co 



2-7 \ ' 



2fcc2fc 2k 
?0 • 



^ 2fc + 2 - 7 



(43) 
(44) 

(45) 

(46) 
(47) 



It is obvious that Q is a positive function of ^o- Therefore, 
from ( p^ one concludes Fx > indicating that M2 is pulled 
away from the centre. M2 will be in static equilibrium if Fx 
is balanced with the gravitational force of Mi, i.e.. 



F, 



GM1M2 GM^ 



(48) 



(2a)2 (2a)2 ' 
By substituting from (^) and (^^ into (^^, we obtain 

S= — ^e~\ s>0, (49) 
roi 

where roi — ro2 , cri = (72 (because we assumed /3 = 0) and 
^^4a^a^-^Q(go)^ (50) 

TrO"! 

Our numerical computations of E reveal that (Too 3> (Ti, 
which guarantees e <^ roi as desired. Following a similar 
procedure as above, one can show that the double nucleus 
remains in equilibrium for /3 > 0. 

5.4 The nuclei of M31 and NGC4486B 

In many respects, the surface density isocontours of our 
model galaxies are similar to the isophotal lines of the nuclei 
of M31 and NGC4486B. Our mass models are cuspy within 
two separatrices. Such curves can be distinguished in the 
nuclei of M31 and NGC4486B (see L96 and T95). We are 
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not sure that the nuclei of M31 and NGC4486B are really 
cuspy, because existing telescopes can not highly resolve the 
regions around PI and P2 (even HST images contain "few" 
bright pixels at the locations of PI and P2). Whatever the 
mass distribution inside these nuclei may be, our models 
reveal that double nuclei can exist even in the absence of 
supermassive BHs. 

The nucleus of M31 can also be explained by the eccen- 
tric disc corresponding to /3 = 1. In such a circumstance, 
stars moving in butterfly orbits form a local group in the 
vicinity of {x = 0,y — 0). The accumulation of stars around 
this local minimum of $ can create a faint nucleus like P2 
(see T95). Therefore, P2 will approximately be located at 
the "centre" of loop orbits while the eccentric, brighter nu- 
cleus (PI) is at the location of the cusp. In other words, 
loop and high-energy butterfly orbits will control the overall 
shape of outer regions, horseshoe orbits will generate PI and 
low-energy butterflies will create P2. 

5.5 Challenging problems 

It is not known for us if there are rotationally supported 
double structures or not. This idea comes from the fact that 
we can replace the point masses of the restricted 3-body 
problem (the restricted 3-body problem is usually expressed 
in a rotating frame) with a continuous distribution of mat- 
ter. Moreover, NGC4486B and the bulge of M31, are three 
dimensional objects and the assumption of planar models 
seems to be a severe constraint. 

So far we showed that the double nucleus can be in 
static equilibrium due to the existing gravitational effects 
of the model. The stability study of such states, however, 
remains as a challenging problem. 

Our next goal is to apply the method of Schwarzschild 
(1979,1993) for the investigation of self-consistency. 
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